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Abstract. Pulse-coupled threshold units serve as paradigmatic models for a wide range of complex systems. When 
the state variable of a unit crosses a threshold, the unit sends a pulse that is received by other units, thereby 
mediating the interactions. At the same time, the state variable of the sending unit is reset. Here we present and 
analyze a class of pulse-coupled oscillators where the reset may be partial only and is mediated by a partial reset 
function. Such a partial reset characterizes intrinsic physical or biophysical features of a unit, e.g., resistive coupling 
between dendrite and soma of compartmental neurons; at the same time the description in terms of a partial reset 
C ~ ) ■ makes possible a rigorous mathematical investigation of the collective network dynamics. 

£N| ' The partial reset acts as a desynchronization mechanism. For N all-to-all pulse-coupled oscillators an increase in 

the strength of the partial reset causes a sequence of desynchronizing bifurcations from the fully synchronous state 
via states with large clusters of synchronized units through states with smaller clusters to complete asynchrony. 

QBy considering inter- and intra-cluster stability we derive sufficient and necessary conditions for the existence and 
stability of cluster states on the partial reset function and on the local dynamics of the oscillators, as specified by their 
i rise function. This analysis of invariance and stability goes beyond monotonicity and curvature of the rise function 

and classifies the local dynamics according to certain expansion and contraction properties that depend also on the 
third derivative of the rise function. We obtain analytical bounds for the stability of cluster states and for a specific 
class of oscillators a rigorous derivation of all N — 1 bifurcation points. Thus the sequence of bifurcation points is 
i extensive (of the order of the system size N) with the actual number of bifurcating states growing combinatorially 

in N . We demonstrate that the entire sequence of bifurcations may occur due to arbitrarily small changes of the 
reset function. We illustrate that the transition is robust against structural perturbations and prevails in presence 
' of heterogeneous network connectivity and rise functions with mixed curvature. 



1. Introduction 

> 

■ Networks of pulse-coupled units serve as paradigmatic models for a wide range of physical and biological systems 

as different as cardiac pacemaker tissue, plate tectonics in earthquakes, chirping crickets, flashing fireflies and 
neurons in the brain [321 EJ HfS U\ ■ In such systems, units interact by sending and receiving pulses at discrete times 
that interrupt the otherwise smooth time evolution. These pulses may be sound signals, electric and electromagnetic 
activations as well as packets of mechanically released stress. Pulses are generated once the state of a unit crosses 
a certain threshold value (e.g. the mechanical stress of a tectonic plate becomes sufficiently large or the voltage 
across a nerve cell membrane becomes sufficiently high); thereafter the state of the sending unit is reset. 

Synchronization of oscillators is one of the most prevalent collective dynamics in pulse-coupled systems [45l [HJ 
E2H[I(2[30l[26l[23[60]. Often not all units are synchronized but form clusters consisting of synchronized sub-groups 
of units which in turn are phase-locked to other clusters [H | [6 ^ [44 t [^[3 ^ [29 t l40 l [50 l [5T]. 

In neuronal networks synchronization and clustering of pulses constitute potential mechanisms for effective feature 
binding. In this paradigm, different information aspects of the same object represented by activity of different nerve 
cells are pooled together by temporal correlations and in particular due to synchronous firing [56], [57]- However, 
strong synchronized firing of nerve cells can also be detrimental: synchrony is prominent during epileptic seizures 
[f8l [42] and observed in the basal ganglia during Parkinson tremor [17]. Here mechanisms for desynchronizing 
neural activity are desirable [59] 141] . 

To study key mechanisms that may underly synchronization, e.g. in biological neural networks, analytical 
tractable models of pulse-coupled oscillator are helpful tools [49] [45] [38l SI 121] [HI ESI HI] ■ Here the rise of the state 
variable of a free oscillatory unit towards the threshold, the unit's rise function characterizes the sub-threshold 
dynamics. If after reception of a pulse the state variable of the unit stays below threshold it is said to receive sub- 
threshold input, whereas excitation above the threshold is supra-threshold. Mirollo and Strogatz [45] showed that 
biological oscillators always synchronize their firing in homogeneous networks with excitatory all-to-all coupling if 
the rise function has a concave shape. The synchronization mechanism they find has two parts: (i) effective decrease 
of phase differences of units due to sub-threshold inputs and (ii) instant synchronization due to supra-threshold 
inputs and subsequent reset to a fixed value. 
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In general, supra-threshold excitation and a subsequent reset is a dominant mechanism for synchronization of 
pulse-coupled oscillators because input pulses that force non-synchronized units to cross threshold at nearby times 
are reset to the same value leaving the units in the same state or in very similar states afterwards. Although this 
reset mechanism plays a crucial role in the synchronization process and the coordination of pulse generation times, 
its implications for the collective network dynamics has - to our knowledge - not been investigated systematically 
so far: most existing model studies reset the units with supra-threshold inputs to a fixed value independent of the 
strength of supra-threshold excitation [6j E0Q GUI HH EU E2 [HI [31] . This results in a complete loss of information 
about the prior state of the units and makes the dynamics non-invertible. Some other studies consider the opposite 
extreme: a complete conservation of supra-threshold inputs during pulse sending and reset [Ml [8]. Here we aim at 
closing this gap by presenting and analyzing a model where the reset (and thus the loss of information about the 
prior state and the strength of supra-threshold excitation) can be varied systematically. 

This article is organized as follows: In section[2]we propose a simple model of pulse-coupled oscillators with partial 
reset, where the response to supra-threshold inputs can be continuously tuned and is not an all-or-none effect. To 
isolate the consequences of partial reset, we focus on homogeneous systems of all-to-all pulse-coupled oscillators with 
convex rise function. In section [3] we briefly present the results of numerical simulations and find that the partial 
reset has a strong influence on the collective network dynamics: Increasing the strength of the partial reset induces 
bifurcations from the fully synchronous state via states with large clusters of synchronized units through states 
with smaller clusters to complete asynchrony. We study this transition rigorously by considering existence and 
stability periodic cluster states with respect to inter-cluster properties in section [5741 and to intra-cluster properties 
in 13.51 We derive conditions on the partial reset function and the rise function to bound regions of existence and 
stability of cluster states. In 13.61 we present a rigorous derivation of N — 1 bifurcation points for a specific class 
of oscillators. We demonstrate that the entire sequence of bifurcations may occur for arbitrarily small changes of 
the reset function, thus underlining the strong impact of partial reset on collective network dynamics. In section 
[4] we numerically illustrate that the transition is robust against structural perturbations and prevails in presence 
of heterogeneous network connectivity and rise functions with mixed curvature. In section [5] we discuss our results 
and relate the simple partial reset model to biophysically more realistic neuron models . 

Specific aspects of the implications of linear partial reset on synchronization properties for oscillators have been 
briefly reported before |35J. 

2. Networks of Pulse-Coupled Units with Partial Reset 

We first propose a class of pulse coupled threshold elements with partial reset and thereafter focus on units that 
oscillate intrinsically. 

2.1. Absorption Rule and Instant Synchronization. We consider TV threshold elements, which at time t are 
characterized by a single real state variable u,i{t) with i G {1,2,..., N}. In the absence of interactions the state 
variables evolve freely according to the differential equation 

(2.1) ^- Ul = F( Ul ) 

with a smooth function F : K — > E specifying the intrinsic dynamics of the units. The free dynamics are endowed 
with an additional nonlinear reset upon reaching a fixed threshold 8 from below 

(2.2) Ui (t~) =6 => Ui{t)=p 

where p < 9 is the reset value and we used the notation m (t ± ) = ]im s \o u i (t± s )- By an appropriate shift and 
rescaling of the state variable and its dynamics we set p — and 9=1 without loss of generality. 

The units are <5-pulse coupled. If unit j reaches the threshold a pulse of strength e\j > is send instantaneously 
to units % and their membrane potential is increased by an amount e,j 

(2.3) Uj (*") =6 => uf ) = Ul (r) + E i:j 

If there is no connection from unit j — > i we set e,j = 0. 
If a unit i crosses the threshold due to a pulse from unit j 

(2.4) Ui (t~) + Eij > e 

it is said to receive supra-threshold input. If a unit receives supra-threshold input it crosses the threshold from 
below, sends a pulse and is reset. Previous models usually reset these units in the same way as if they reached the 
threshold without this recurrent input. In this type of reset, also referred to as the absorption rule (e.g. [45] ) 



(2.5) 



Ui (t)>9 => u t (t+) = p 



a b 




Figure 2.1. Model dynamics, (a) Sample traces of three units with (i) network connectivity Sji = 
Ski > 0. (ii) At time t = t\ unit i reaches the threshold 9 and its membrane potential is reset to p. It 
generates a pulse which is send to the units j and k. (iii) unit j receives the pulse and its membrane 
potential is increased to Uj(t^)+£ji, the pulse is sub-threshold, (iv) unit k receives a supra-threshold 
pulse, Uk{t\) + £kj > 0, its membrane potential is set to R(() = R (itj(ij~) + £kj — 9) using the 
partial reset function R. (b) Sample avalanche process with 6 = {1, 2, 3} and n = 3 in (i) a N = 4 
all-to-all network £y = (1 — 5y)e. (ii) unit i = 1 reaches the threshold O^ ' = {1} and sends a 
pulse to the other units (arrows). Their potentials are updated to = Ui(t~) +£n, causing unit 
i = 2 to generate a pulse, &^ = {2}. (iii) The pulse is received by the other units yielding a 
potential = u'p + ea which brings unit i — 3 above threshold O^ 2 ^ = {3}. (iv) The potentials 
become u\ = uf^ + e& and no further unit crosses the threshold, = 0. (v) The avalanche 
stops and units that received supra-threshold input are reset to Ui(t) = p + R (itf — 9 J . 



the total supra.threshold input is lost. As a consequence two or more units initially in different states m and 
simultaneously receiving supra-threshold inputs will all be reset to the same value p, making the absorption rule 
a strong instant synchronizing element of the network dynamics. An alternative considered in previous studies as 
[38| [31] is total input conservation, 

(2.6) m (t)>9 Ui (t+) =p+( Ui (t) - 6>) 

i.e. the total supra-threshold input charge C = Ui{t) — 9 is added to the potential p after the reset. 

2.2. Partial Reset. Here we propose a more general model where the reset value is given by a partial reset function 
R(C) that depends on the supra-threshold input charge C = Ui(t) — 9, 

(2.7) Ui(t)>6 Ui(t + ) = p + R(ui(t) - 9) 

We assume that supra-threshold inputs only lead to excitatory contributions and thus define: 

Definition 1 . A function R : R — * R which is monotonically increasing and satisfies R (0) = is called a partial 
reset function. 

For a linear partial response we set 

(2.8) R c (C) = c( 

with the remaining fraction < c < 1 of supra-threshold input charge after the reset. For c = we recover the 
absorption rule (|2.5p while c = 1 corresponds to total charge conservation l|2.6p . 

Motivation for this extension comes from neural networks. Neurons consist of functionally different compart- 
ments, including the dendrite and the soma. While synaptic input currents are collected in the dendrite, the 
electrical pulses are generated at the soma. Additional charges not used to excite a spike may stay on the dendrite 
and contribute to the membrane potential after reset at the soma. Due to intra-neuronal interactions a part of this 
supra-threshold input charge may be lost due to the reset at the soma: 



Definition 2. A partial reset function R is said to be neuronal, if < R(() < C f° r a U C > 0- 

2.3. The Avalanche Process. Since the interaction is instantaneous, a pulse generated by unit j may lift other 
units above threshold simultaneously. These then generate a pulse on their own, etc. This leads to an avalanche of 
pulses (cf. fig. 12. ip : Units reaching the threshold at time t due to the free time evolution define the triggering set 

(2.9) 0(°) = {j | Uj (r) = 0} 

The units j S O^ generate spikes which are instantaneously received by all the connected units i in the network. 
In response, their potentials are updated according to 

(2.10) uj 1} := u % (r) + J2 su 

iee 

The initial pulse may trigger certain other units k 6 O^ 1 ' = |fe | Uk (t~) < 1 < to spike, etc. This process 

continues n < N steps until no new unit crosses the threshold. At each step m S {2, 3, . . . , n} the potentials are 
updated according to 

(2.11) := + J2 *« 

jee m 

where 

0< TO > = {k\u { ™- 1] < 1 <4 m) } 
The potentials immediately after the avalanche O = Ug=o °^ s ^ ze a ~ I® I are obtained via 



(2.12) u % (t+) 



'ui(t )+J2 je e £ v 1 £ 9 



using the partial reset R for units having received supra-threshold inputs. 

In general there is an ambiguity in fixing the precise order of potential updates and resets during the avalanche. 
For example an instantaneous reset could directly follow a single supra-threshold input pulse instead of resetting the 
potentials at the end of the avalanche. Motivation for our choice again comes form neuroscience: It is valid when 
the time scale of the action potential (and subsequent reset) is much faster than the time scale of the synaptic input 
currents. These in turn should be much faster than the time scale of the mechanism reducing the supra-threshold 
input, e.g. the refractory period (cf. the forthcoming publication |36j). Our model (|2.12[) then is the limit where 
all these time scales become small compared to the time scale of the intrinsic interaction free dynamics. 

For non-zero partial reset functions potential differences of oscillators involved in a single avalanche will in general 
not be fully synchronized after the reset. Thus despite the fact that units are generating pulses simultaneously they 
can have different phases afterwards. We therefore distinguish between phase synchrony where units have identical 
phases and the weaker condition of pulse synchrony which corresponds to simultaneous firing only but allows 
differences in the phases. When examining the system with a higher time resolution phase synchronized units will 
stay synchronized whereas pulse synchronized units fire within a short time interval. 

2.4. Phase Representation of Pulse-Coupled Oscillators with Partial Reset. In the remainder of this 
article we will concentrate on units with strictly positive F > in (|2.1j) . Then the individual units become 
oscillatory as the strictly monotonically increasing trajectory Ui(t) of a unit i starting at Ui(0) = reaches the 
threshold after a time T and is reset to zero again. By an appropriate rescaling of time we set T = 1. Defining a 
phase like coordinate (cf. [45]) via 

/•««(*) i 

(2.13) & (t) = U- 1 ( Ui (t)) := / —— -du 

Jo F ( u ) 

the interaction free dynamics simplify to 

(2-14) j t Ut) = 1 

By definition U~ x is strictly monotonically increasing and has a strictly monotonically increasing inverse U. By 
our choice of normalization they obey XJ~ X (0) = = U (0) and U^ 1 (1) = 1 = U (1). Note that the function U ((f)) 
captures the intrinsic rise of the membrane potential towards the threshold. 



Definition 3. A smooth function U : [0, oo) — > [0, oo) is called a rise function if it is strictly monotonic increasing 
U' > and is normalized to U (0) = and U (1) = 1. 



Definition 4. Given a rise function U and a partial reset function R we define for e > the (sub-threshold) 
interaction function H £ : [0,oo) — > [0,oo) by 

(2.15) H e (0) := H (</>, e) := U' 1 (U (0) + e) 
and the supra-threshold interaction function J e : [U^ 1 (6 — e) , oo) — * [0, oo) by 

(2.16) J e ((f)) := J ((f), e) := U' 1 (R (U (</>)+£- 9)) 

We remark that H^ 1 — H— e , The pulse-coupling in the potential representation eq. Q2.3jl carries over to the 
phase picture using the interaction function H: 

(2.17) &(t~)=l => <j> i (t + T)=H{4> i ({t + T)-),e ij ) 
Equation l|2.12jl for phases after an avalanche 8 at time t becomes 

(2.18) , ((t+ )H ff K' ,E ^ '** 

3. Network Dynamics 

To identify the effects of the partial reset on the collective network dynamics we here first focus on homogeneous 
networks consisting of TV units with all-to-all coupling and without self-interaction, i.e. 

(3.1) e lJ = {l-8 r] )e 

i,j € {1, 2, ... , N}. We impose the condition J^j £ ij = (N ~ 1) £ < # — p = 1 to avoid self-sustained avalanches of 
infinite size. 

The analysis of Mirollo and Strogatz [45J shows that in this situation (with a slightly different avalanche process) 
synchronization from almost all initial conditions is achieved when the rise function is concave (U" < 0) and the 
absorption rule (R = 0) is used. In fact, their result can be generalized to the partial reset model used here and any 
partial reset function R that is non-expansive (e.g. R' < 1). For expansive R the synchronized state no longer has to 
be the global attractor of the dynamics and typically irregular dynamics are observed. The proof of synchronization 
for non-expansive partial resets and the analysis of the bifurcation to the irregular dynamics are presented in a 
forth coming study [36] . 

In this article we concentrate on convex rise functions U, i.e. 

(3.2) W U(<P)> °- 

This property holds for large class of conductance based leaky-integrate-and-fire (LIF) neurons and a class of 
quadratic- integrate-and-fire (QIF) neurons (cf. appendix (jBj ) . Studying convex rise functions is further motivated 
by the fact that for these rise functions we already observe a rich diversity of collective network dynamics with a 
strong dependence on the partial reset R. However, some of our results also apply to more general rise functions 
and in particular to sigmoidal shapes as often found for neurons [191123] . 



3.1. Numerical Results: A Sequence of Desynchronizing Bifurcations. Systematic numerical investigations 
indicate a strong dependence of the network dynamics on the partial reset R. In particular, we find synchronous 
states, cluster states, asynchronous states and a sequential desynchronization of clusters when increasing the partial 
reset strength, e.g. by increasing the parameter c when using R = R c . 

Starting in the synchronized state and then applying a small perturbation to the phases we observe that the 
synchronized state is stable for sufficiently small c (fig. 13.11 1). When increasing the partial reset strength the 
synchronized state becomes unstable and we observe smaller clusters in the asymptotic network dynamics (fig. 13.11 
II) where the final cluster state depends on the precise form of the perturbation. The maximally observed cluster 
sizes depend on the value of c (fig. I3.3|l . For sufficiently large c only the asynchronous splay state, i.e. a state with 
maximal cluster size a = 1 , is observed (fig. 13.11 III) - 

Starting from random initial conditions we find that for sufficiently small c the synchronous state coexists with a 
variety of cluster states and the asynchronous state (fig. 13.21 1 and fig. I3.3[l . Increasing c, the states involving larger 
clusters become unstable (fig. 13.21 II) until finally all random initial conditions lead to the asynchronous state (fig. 

El in). 

What is the origin of this rich repertoire of dynamics and which mechanisms control the observed transition 
of sequential desynchronizing bifurcations? To answer these questions, we analytically investigate the existence 
and stability of periodic states involving clusters of arbitrary sizes a < N. The following analysis reveals that 




Figure 3.1. Desynchronization transition in a network with parameters N = 50, e = 0.0175, 
U = Ub, b = —3, R = R c for different values of the partial reset strength: (I) c = 0.025 < c C i , 
(II) c = 0.5 € ^Ccf' 1 , Cc?'^ and (III) c = 0.7 > Cc?. Plotted are (a) the phases fa (t r ) of all units at 

pulse generation times t r of the r-th spike of a reference unit (cf. return map (GLU), (b) potential 
traces of all units and (c) the raster plots marking the times of pulse generation of each unit i. 
The network is started in the synchronous state and then a small perturbation is applied at a time 
indicated by arrows. 



the sequence of bifurcations is controlled by two effects: sub-threshold inputs that are always synchronizing and 
supra-threshold inputs that are either synchronizing or desynchronizing depending on the partial reset strength. 

3.2. Strategy of the Analysis. We split up our analysis of the dynamics into two parts. First we assume that 
all avalanches are invariant, i.e. the given clusters do not decay into smaller sub-clusters. This assumption allows 
to group all oscillators firing in a single avalanche together into a single "meta-oscillator" with increased firing 
strength and an effective self interaction (cf . fig. I3-4|) . The analysis of a state $ in homogeneous all-to-all network 
(eij = (1 — Sij)e) of TV oscillators with avalanche sizes a s , s£ {1,2,..., to}, J2 s a s — N then reduces to analyzing 
a network of to meta-oscillators with coupling strengths 

(3.3) Eij = (1 - Sij)ei + SijSu 

and Si = a,i£, en = (a t - l)e. 

In a second step we derive conditions under which an avalanche of a certain size will indeed not decay into smaller 
groups. 

3.3. Notations: State Space, Firing and Return Map. A state of a network of N pulse-coupled oscillators is 
completely specified by a phase vector 

(3.4) <f> = (fa,fa,...,<f> N )eS N = S 1 x- - -xS\ 

N times 
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Figure 3.2. Same as in Fig. I3.1l but starting the dynamics form random initial phases distributed 
uniformly in the interval [0, 1). (Illb) As the transient is similar to (lib) a close up of the asymp- 
totic asynchronous (splay) state is shown instead. In (c) the oscillator are labeled by increasing 
initial phases. 



where 0; S S 1 = R/Z are the phases of the individual units. Since the time evolution in between avalanches is a 
pure phase shift (|2 . 14[) it is convenient to consider a Poincare section S of S N with states just before the firing of 
one or more oscillators, i.e. 

(3.5) S={<PeS N \33e{l,...,N}, ^- = 1} . 

It is convenient to relabel the oscillators after each avalanche such that 1 = <j>\ > $i > ■ ■ ■ <j>N-i > <Pn > 0. To 
specify the state of the network completely the permutation 7r _1 used for relabeling of the oscillators is remembered. 
The largest phase tfii = 1 thus belongs to the oscillator i = 7r(l), the second largest 4>2 to i = 7r(2), etc., and the 
smallest (f>N to oscillator i — tt(N). Thus an equivalent description of the state space S is given by 

(3.6) S p = {((0 2 , . . . , <f> N ) , tt) € S N - X X S N \l>(f> 2 > ■■■ <f>N-l ><f>N>0} 

Here Sn is the group of all permutations of N elements. We use the convention that all index labels i are taken 
modulo the network size N, e.g. labels i = j and i = N + j denote the same oscillator. 

The Poincare map of the network dynamics for the Poincare section S is the firing-map K that maps the state 
$ € S of the network just before the s-th firing time t s of an avalanche to the state just before the next avalanche 
that occurs at time t a +i* 

(3.7) K ($ (t-)) = $ (t- +1 ) € S 

Having determined the next avalanche & from a state $ £ S, the map K is a composition of the avalanche map 
(|2.18p and a subsequent shift of all phases to a state in S. Note that the firing map is fully determined by the pair 
(Q, a) which is a function of We denote a phase shift of size a by 




Figure 3.3. Sequential desynchronization transition in a network of N — 50 units (U = Ub, 
b = —3, e = 0.0175). Asymptotic cluster sizes a observed in the asymptotic network dynamics of 
500 simulations for each c G {0, 0.025, . . . , 1} starting from a perturbed synchronous state (•, cf. fig. 
I3.ip or from other random phases distributed uniformly in [0, 1) (o cf. fig. 13 . 2[) . The shaded area 
marks the sequential desynchronizing transition, the solid line shows the exact theoretical prediction 
(|3.36|) continuously interpolating the N — 1 bifurcation points cit\ a £ {2, 3, ... , N} above which 

(2) 

clusters of size a become unstable. The gap for clusters sizes 43 < a < 50 at < c < Ccr appears 
as cluster states involving these avalanche sizes do not exist according to lemma [71 




Figure 3.4. Network reduction method, (a) A network of N = 4 oscillators producing an 
avalanche of three units as in Fig. 12.11 may be effectively reduced to (b) a smaller network of 
size N = 2 with inhomogeneous coupling and self-interactions by grouping the oscillators involved 
in the avalanche to a single meta-oscillator. The effective network connectivity then has the form 
lEOl) . 



(3.8) S (0, a) := S a (0) := <\> + a 

The equivalent firing map acting on the state space S p is denoted by K p . For the phase part we write 

To track the network dynamics we consider a mapping of the state just before a fixed reference oscillator k fires in 
an avalanche at time t r to the state just before this oscillator fires again at t r+1 : 

(3.9) M ($ (t-)) = $ (t- +1 ) 

M is called the return-map and is the Poincare map of the system on the section {$G6>|0fc = l}. Again the 
equivalent return map acting on S p is denoted by M p . The number m of avalanches occurring in the application 
of the return map is a function of the initial phase vector $ = $ (t~) and the return map M then is a composition 
of m firing maps K. Thus M is completely specified by an ordered firing sequence 

(3.10) ^ = ^($) = {(e s , CTs )}r =1 

where the pairs (6 s ,er s ) specify the avalanche set 6 S and subsequent shift a s of the s-th firing map. 



Given a firing sequence p.lOp . we set a s = \Q S \ and in the case of homogenous networks with coupling (|3.ip . 
e s = a s e. A composition of shift and interaction maps is denoted as 

m 

(3-11) O (S* s o H £s ) (fa := S am o H Sm o S 0m _ t o H^... o S CT2 o H £2 o o ff El (0) 

8=1 

3.4. Existence and Stability of Asynchronous Periodic States in Meta-Oscillator Networks. 

Definition 5. An asynchronous periodic state of a network of N pulse-coupled oscillators is a state $ G S which 
is invariant under the return map, i.e. M ($) = $, and with avalanche sizes a s = 1, s e {1,2,..., A 7 "}, i.e. each 
oscillator generates a pulse separately and does not receive supra-threshold input. 

Initially assume that all clusters stay forward invariant, i.e. do not decay in to smaller sub-clusters during the 
network dynamics (cf. sec. I3.2|) and thus consider networks of meta-oscillators with effective coupling matrix (|3.3|) . 
A periodic cluster state in the original model thus becomes an periodic asynchronous state in the reduced effective 
meta-network. In the following we derive conditions for the existence of the asynchronous state and its stability in 
a meta-network. 

Lemma 6. Consider a network (|2.17|l - i|2.18p of N oscillators with pulse coupling matrix (|3.3p and neuronal partial 
reset. Let S = (oi, . . . , a N ) G R N and define L:R N x S 1 by 

N+i—l 

Li(E,0):= O (S<r.oH e .)oS n oJ Bti (<l>) 

s=i+l 

for i £ {1,2,..., N}. Then the asynchronous state exists if and only if there is a solution E* G ]& N to the equation 
(3.12) L(S,1) = (1,1,..., 1) 

that satisfies a* > for all r G {1, 2, . . . , N}. 
Proof. Assume there is a solution £*, a* > 0. Set 

K = i, <t>l = (h- 1 ° s-}) o (jj- 1 o s-}) o . . . o (h-\ o s-}_) (i) 

for i G {2 . . . , N}. Using that E* is a solution to (|3. 12|) we have (S a * o £T er ) o S^* o J eNN (l) = 1 and 

0tv = S a ' o J £JVJV (1) > since a* N > 0. Further using Ei > and er* > the phases are ordered according to 
<t>*x = 1 > <t>*2 > ■ ■ ■ > 4>* N > and $* = . . . , 4>*n) € 5. 

Starting from the state <&* the first pulse of oscillator i = 1 results in potentials = i? (£11), u| = U ((f)*) +£i, 
i € {2,3, . . . , A 7 }. Since i2(e«) < £« <£,-<£» + C/(0) for all > and H Ei (fa < H Ei (tp) for < ^ we have 
< < < • ■ ■ < Further 

4 X) = ^ (#) + ei = ((h- 1 o (1)) + £i = J7 (1 — ^) < 1 

as cr* > 0. Thus oscillator i = 1 fires without triggering any further oscillators yielding an avalanche set 9i = {1}. 
In addition the oscillators have to be shifted by to obtain fa = 1. Thus the first pair in the firing sequence is 
({l},er*). Applying the same arguments to the new phases $*W = K(<f>*) yields ({2} ,a*) for the second pair. 
Repeating these steps N times one obtains a firing sequence 

^($*) = {(W,<)}f =1 

Thus 

N i-1 

Mi(**) = O (S (T ,oH er )oS at J eu oQ(S K oH er )(fa) 

r— i+l r—1 
N 

r=i+l 

= (H- 1 o S-}) o (if-* o ^) o • • • o (ff-^ o g-lj (1) 

= 01 

using p,12p in the third row. Hence M ($*) = $* and the asynchronous state is invariant under the return 
map. 



Conversely an periodic asynchronous state yields a solution to (|3. 12|1 . since by definition no oscillator receives 
supra-threshold input and thus there are a phase shifts <7j > after each pulse generation of oscillators i <E 
{1, 2, . . . , N}. Invariance of the periodic asynchronous state then shows that in fact £ = (ctj., 02, . . . , &n) is a 
solution to (J3TT2J) . Hence there is no periodic asynchronous state if the solution does not exist. If there is a solution 
with a* < 0, let s be the smallest index such that er* < 0. Starting in the state $* the first firing of oscillator i = s 
will cause oscillator i = s + 1 to fire in the same avalanche since its potential at this point is given by 



«& = U 



r=l 

U [H- 1 o S~} (1)] + s s = U(1- <j* s ) > 1 



i.e. {r, r + 1} C Q r and the system is not in a periodic asynchronous state. □ 

Corollary 7. In a network (|2 . 1 T[) - (|2 . 18[> of N oscillators with homogeneous all-to-all coupling matrix (|3.1|l an 
asynchronous (splay) state exists. 

Proof. Let 

N-l 

L(a) := O {SaoH e )oS a oJ (l) 
s=l 

Now since L(0) = U^ 1 (e(N — 1)) < 1 and ^L(cr) > 1 the intermediate value theorem ensures the existence of a 
a* > satisfying L(a*) = 1. S* = (a*, . . . , a*) is a solution to p,12p . This result is independent of the partial 
reset R as en — and no oscillator receives supra-threshold input in the asynchronous state. □ 

In Fig. 13.31 we observe no cluster states involving avalanches of size 43 to 49. This is precisely because f3. 121 has 
no solutions when setting £j = a^e, en = (a^ — l)e for a± € {43,44, . . . ,49} and any further < a, G N, i > 1 and 
m such that fls = 

Note that lemma [7] holds for any rise function U. If there are q different positive solutions to (|3. 12|) there coexist 
q different periodic asynchronous states. A convex U ensures that the solution is unique because L (E, 1) then 
becomes invertible for all S 6 R N . Another consequence of convexity is that given the existence of an asynchronous 
state in a meta-oscillator network it is linearly stable as the following theorem shows: 

Theorem 8. Consider a network (|2 . 1 T|) - <|2 . 18|1 of N oscillators with pulse coupling matrix l|3.3p and neuronal 
partial reset. If a periodic asynchronous state exist it is linearly stable. 

Proof. Existence of the asynchronous state ($*, id) € S p with <I>* = (0| 5 • ■ ■ , 4>*n) implies invariance under the return 
map M p , 

(3.13) M p ($*,id) = (**,id) 

For the intermediate states we set 

($W,7rOO) := (K?) s ($*,id) se {0,1,2,.... N} 

If oscillator i generates a pulse all oscillators j ^ i receive the same input Ej and oscillator i receives an input 
Ejj < £j, Hence, using R(() < (, we find that the oscillators do not change their firing order and it^ is a cyclic 
permutation to the left 7P S ) (i) = i — s. 

We show that the asynchronous state is linearly stable: Adding a perturbation A( ' = (s^, ^-1) *° ^ ne 
asynchronous state such that initially the phases are given by 

*W := (^°\... J 0W 1 )=** + AflO 

We take the perturbation to be sufficiently small such that the oscillators still fire asynchronously, i.e. the avalanches 
are of size a s = 1 and the order of the events is preserved. In the following, terms of order O ^(A^ ^) 2 ^ are neglected 
which we indicate by a dot above the equality sign (=). After s firing events the phases are 

where 



is the phase perturbation before the next firing and A^) is the Jacobian matrix of K|, at ($^ s 1 \7r^ s l n 



1 — H (-02, £ tt(i)) the phase part of the firing-map for N > 3 is 



(3.14) 
Setting a 

(3.15) 

Inserting (|3.15p into (|3. 14|) gives 

(3.16) 

with 



.4 



is) 



d(f)j 



L ($0-1) tt^-I) 



/ H (^3,^(1)) + C \ 

H (04,^(1)) + CT 
\ J (l,£7r(l)7r(l)) +0" / 



AW = 



l 2 °3 




(») 

2 

(s) 
2 



o \ 



l (s) 

N 



(3.17) 



u 



^'K (<^~ 13 )) 

Since £j > it follows that H E] (0) = L^ 1 (U(<p) + £,-) > <f>. Thus a| s) < 1 since J7 is convex. Also U' > and 



hence 
(3.18) 



< a^ s) < 1 



Now the Enestrom-Kakeya theorem (cf. appendix [A] and [24l [341 [4j [32]) applied to the matrix A^) shows that with 
these properties the spectral radius p (A' s )) of A' s ) satisfies 



max 

ie{l,—,N-i} 



a« < 1 



Thus 
(3.19) 



A (nN) 



N 



\r=l 



0) A (0) 



AT 
r=l 



as n 



and the asynchronous state is linearly stable. For N = 2, p (A^)) = a% < 1. 



□ 

This result is illustrated in Fig. 13.51 Due to the convexity of the rise function oscillators perturbed to larger 
(smaller) phases compared to the asynchronous state are less (more) advanced by input pulses pulling the perturbed 
phases back to the invariant asynchronous dynamics. 

Combining corollary [7] and theorem [8] we obtain: 



Corollary 9. In a network (|2 . 1 7[) - (|2 . 18[1 of N oscillators with homogeneous all-to-all coupling matrix Q3.1|l . neu- 
ronal partial reset R and convex rise function U the periodic asynchronous (splay) state exists and is linearly stable. 

3.5. Impact of Partial Reset on Intra-Cluster Stability. In the state of synchronous firing all units in an 
all-to-all coupled network receive a supra-threshold input pulse of strength (N — 1) e suggesting a rather strong 
influence of the partial reset R onto the network dynamics. Indeed, as shown in Fig. 13.31 for the partial reset R c 
one observes a sequential destabilization of clusters when increasing the reset strength c. In this subsection we 
study this behavior analytically and explain the observed transition. The strategy is to focus on a single cluster 
and derive general conditions which ensure the stability of this cluster under the return map. As the return map 
depends on the firing sequence T only we give lower bounds below which cluster sizes are ensured to be stable and 
upper bounds above which cluster sizes are unstable. However, we find that for a special class of rise functions a 
full analytical treatment is possible. 




Figure 3.5. Stability of the asynchronous state, (a) Graph a network of N = 2 oscillators 
with connectivity £y = (1 — 6ij)£j with < £2 < £i- (b) Firing of oscillator i = 1. For the 
oscillator i — 2 with initial phase — <f>^ + 8^ smaller than in the invariant asynchronous 
state (j)^ (gray) the input advances the phase ip^ more in comparison with the advance of (f)^ 
in the asynchronous state due to the convexity of the rise function, (c) After the interaction a 
subsequent shift completes the firing map K. In total the derivation from the asynchronous state 
8^ has become smaller, (d) Firing of oscillator i = 2. Phases which are perturbed to larger values 
than the asynchronous state are less advanced by inputs due to convexity of the rise function, 

i.e. H £l ^V , 2°' > ) — V4°^ > Hei yf^J — 4 > 2 > ' > ■ (A) I n total the return map M decreases the phase 
perturbations \S^\ < \S^\. These stabilizing dynamics of the asynchronous state due to the 
convexity of the rise function generalizes to larger networks as proven in theorem [HI 



Definition 10. A firing sequence T is admissible if there is a state $65 which has firing sequence T = 
It is further called trigger invariant if for the oscillators i € 6 X = {j G {1, 2, . . . , N} | <pj = 1} triggering the first 
avalanche of the state $ = (0i, . . . , 4>n) (cf. I|2.9p ) the return map satisfies Mj($) = 1. Thus for a trigger invariant 
firing sequence T with m intermediate avalanches <df ^ C 8^ +1 . The set of all trigger invariant firing sequences is 
denoted by T. The set of T € T with initial avalanche size ai = is denoted by T ai . 

Let us focus on a single avalanche of size a\ in the network dynamics. To ensure that all units in this avalanche 
fire together again after the return map is applied all units in this avalanche which were triggered to fire by 
a € {1, 2, . . . , a\ — 1} preceding spikes i.e. with phases in 

7j = [U- 1 (1 - as) , 1] 

have to be triggered again after applying the return map. Given a firing sequence T = {(e r , a r )}™=i the return 
map for oscillators i € Oi in the first avalanche is given by 

m 

M r {4>) ■= O ( S °r ° H *r) ° S ai O J ei (0) 
r=1 

Hence the conditions 

(3.20) Mr (4 T ) C II 

for all a S {1, . . . , a% — 1} and all admissible firing sequences T S T ai ensure a cluster of size a% to not split up 
under return. By finding the most synchronizing and most desynchronizing firing sequences (i.e. "best" and "worst 
case" firing sequences) T £ T ai these conditions yield upper and lower bounds for stability of a cluster of size a\ 
under the return map. 



Lemma 11. Consider a network i|2.17p - (|2.18p of N oscillators with homogeneous all-to-all coupling matrix l|3.ip . 



Set 

K 1 = ^ M r (U- 1 (1 - as)) 

and 

bl 1 = sup M r (Z7 _1 (1 - as)) 

J r £T ai 

Then the conditions 

(3.21) w a a x > U' 1 (1 - as) 
for a € {1, 2, . . . , a\ — 1} are sufficient and 

(3.22) ft* 1 > C/" 1 (1 - as) 

are necessary for a cluster of size ai to be invariant under return. 

Proof. ^Mjr(^) > and thus conditions (j3.20|) are equivalent to Mr (U^ 1 (1 - as)) > U^ 1 (1 - as) for a S 
{1, 2, . . . , a\ — 1} and all admissible J 7 E T ai . □ 

Finding the wj 1 and for general U and i? can be done numerically using optimization techniques. However, 
there are two classes of rise functions which allow further analytical investigation. Most of the commonly used rise 
functions, as e.g. the rise function of the LIF neuron or the conductance based LIF neuron fall into one of these 
classes (cf. Appendix [Bj . 

Two oscillators initially at phases <fi and (j> + A(j> receiving a pulse of strength s will have a new phase difference 

(3.23) AH {(/), A<t>, e) := H E (0 + A0) - H s (</>) . 
We denote the domain of AH as 

V := {((t>,A(j),s) |0 < s < 1, < <p < 1, < Ac/) < U' 1 (1 - s) - <j)} . 
Definition 12. A rise function U is increasing the change of phase differences (icpd) iff 

d 

(3.24) — AH U, A<f>, s) > for all U, A<j>, s) € V . 

dtp 

Conversely, it is decreasing the change of phase differences (dcpd) iff 

d 

(3.25) — AH U, A<p, e) < for all (<j>, A<f>, s) € V . 

Off) 

As shown in appendix IB .21 the icpd (dcpd) property is related to the third derivative of U. 

The following lemma allows to bound the change in phase differences if the rise function is icpd or dcpd: 
Lemma 13. Let s r , a r > 0, r € {1, 2, . . . , m}, e — J2T=i £ r> a i > 0. Choose a <j u > such that 

m 

0(S ar oH er )(<f>)<H e oS au (<f>) . 

r=l 

and let tp < cf>. Then for an icpd rise function U 

m m 

(3.26) S ai o H E (0) - S ai o H E (ip) < O (Sa r ° H £r ) (0) - Q {S ar o H Er ) < H E o S au {4>) - H E o S« u (V>) 

r— 1 r— 1 

// U is dcpd then 

m m 

(3.27) S ai o H E (</>) - S ai o H E {i>) > O (S* r ° # Er ) (0) - Q ( S <* ° ff »r) WO > H e ° S* v (0) - H E o S ffu (V) 



Proof. Consider icpd rise functions first: To show the first inequality of eq. (|3.26|) we use induction on m. The 
statement is clearly true for m = 1. Assume it is true for m > 1 then 

S ai o H £ (0) - S^oH, (VO = H £m+1 o H £ _ £m+1 (<t>) - H Sm+1 o ff £ _ £m+1 

= AH (ff e - £m+1 (V') , #e-e m+1 (0) - ffe- £m + 1 (V) , £m+l) 



<AH H £ _ £m+1 , O o ff £r ) (0) - Q o ff £r ) ty) , £m+1 

\ r=l r=l / 

Cm m m ^ 

© {S ar o ff £r ) (V.) , O (S CT „ o H Sr ) (0) - Q (5^ o H Er ) (</>) , £m+l 

r— 1 r— 1 r— 1 y 
m+1 m+1 

= O ( S »> ° ^r) fa) " O ° W 



r— 1 r— 1 

where we used the induction hypothesis and -t^AH > (cf. Ij3.23p ) in the third, and in the fourth again the icpd 

property and the fact that H e _ £m+1 (i/j) < Q" l =1 (S ar o H er ) if YlT=i £ r = e i a i — 0- Substituting < with > we 
obtain the result for dcpd rise functions. 

For the second inequality we also use induction over m. The statement is trivially true for m = 1. Let it be true 
for m > 1 and let a u > such that ©"lY (5 CTr o fl" £r ) (0) < i? £ o S a „ ((/>). Then 

h £ o s au {</>) - # £ o s au y>) = H £m+1 o # £ _ £m+1 o s au (0) - H £m+l o # £ _ £m+1 o s ff „ 

= Ai? (ff e _ £m+1 o fofr) , ff e - Em+1 ° S au (0) - F e _ Em+1 o S au (V) , e m+ i) 

(mm \ 
H e -e m+1 O (V) , O ( 5 <v ° W ~ O ^ ° H *r) W . £ m+l 

r=l r=l / 

(m m m \ 

O (Sa r ° ff £r ) (VO , O (S CTr o H Sr ) (cf,) - O ° H *r) W > Zm+l 

r—1 r— 1 r— 1 / 

m+1 m+1 

= O ° H *r) (<t>) - O ^ ° M 

r—1 r—1 

where in the third row we used the implication 

m+1 m 

(3-28) O ^ ° H e° ) fa) ^ ^ ° fa) O ( 5 *. ° ) M ^ ° ^„ fa) 

s=l s=l 

to apply the induction hypothesis. In the fourth row we again used g^AH > 0, eq. (|3.28p and the icpd property. 
Substituting > with < we obtain the result for dcpd rise functions. □ 



This result is illustrated in Fig. 13.61 (a-c) for icpd rise functions. 

The next theorem allows to determine bounds on the network parameters in order to ensure invariance or decay 
of avalanches under return: 

Theorem 14. Consider a homogenous excitatory all-to-all network of N pulse- coupled oscillators evolving according 
to i|2.17p - l|2.18p with neuronal partial reset R. 
For icdp rise functions U the conditions 

(3.29) I!' 1 (R {(ax - l)e)) - U' 1 (R ((a x - l)e - as)) < U' 1 (1-(N- ai)e) - U' 1 (1 - (N - o x )e - as) 

for all a G {1, 2, . . . , ai — 1} are sufficient to ensure the invariance of an a\-avalanche under return. Necessary 
conditions are 



(3.30) If- 1 (R ((ai - l)e) + [N - ax)e) - U' 1 (R ((ai - l)e - as) + (N - ai )s) < 1 - U' 1 (1 - as) . 

Likewise for dcpd rise functions U sufficient conditions are (|3.3Q[) and necessary conditions l|3.29p for an a\- 
avalanche to not split up under return. 

Proof. Using lemma [13] we find for an icpd rise function and T € T ai 




o io 4> io <i> 1 

Figure 3.6. Rise functions with increasing change (icpd) and no change of phase differences, (a-c) 
Icpd rise function. An initial phase difference changes to A^ after applying a combination of 
interaction maps H £r (blue) of total strength e — J2"=i £r anc ^ shifts S ar (green) such that the final 
maximal phase values are identical, (a) For icpd rise functions the difference A^ is the smallest 
when the interaction is applied in total before the shifts, i.e. H e o S ai and (c) largest when applied 
after the shifts S au o H e . (b) All other maps ©™ =1 (H £s o 5<T a ) produce phase differences which lie 
in between these extremal values (cf. lemma [T3j) . (d-f) The rise function Ub is icpd and dcpd, i.e. 
the phase difference A(f>^ is independent of the order in which the interactions and shifts are 
applied. 



Af>(l) - Mr (CT 1 (1 - as)) = Q {S a „ o H £r ) (S ai ° J £l (1)) - Q (S ar o H £r ) (S ai ° J £l (U' 1 (1 - oe))) 

r=2 r=2 

< H (N _ ai)£ o S„ u o J £l (1) - H (N _ ai)£ o S„ u o J £l (p- 1 (1 - as)) 
= 1 - H (N _ ai)£ o S a o J £l (U- 1 (1 - as)) 

with 

(3.31) a u = U- 1 (1 - (N - oi) e) - U~ l (R (ai - 1) e) 
and thus 

= H {N _ ax)e o o J £l (U- 1 (1 - oe)) 
in l|3.2ip yielding conditions (|3.29|) . Similarly we find 

= 1 - H (J v_ ol)E o J £1 (1) + i? (A r_ ai)e o J £1 (U- 1 (1 - ae)) 

which yields the necessary conditions (|3.30|1 . For dcdp rise functions the expressions for w^ 1 and b^ 1 are interchanged. 

□ 

Proposition 15. In a homogeneous all-to-all coupled network of N neural oscillators evolving according to (|2 . 1 T[) - 
(|2.18p with neuronal partial reset R and dcdp rise function U the condition 

(3.32) R'(C) > ui{u Z^_\n-iI + ) for al1 - ^ ^ ~ ^ 
is sufficient to ensure non-invariance of an a\-cluster under return. 

Proof. Using lemma[13]for a given firing sequence T = {(a s e, cr s )}™ =1 we estimate 

(3.33) 1 - H {N _ ai)£ o S au o J ei (1 - Acj>) < 1 - Mjr (1 - A</>) < 1 - S ai o H {N _ ai)e o J £l (1 - A</>) 

with a u as in l|3.3ip and 07 = 1 — i?(jv-ai)e (^aie(l))- In general a ai-cluster is triggered by a single oscillator and 
thus if 

(3.34) 1 - H {N _ ai)£ o S au o J ai (1 - Ac/)) > A<f> 



for all < A<fi < 1 — U^ 1 (1 — e), we have 1 — Mjr (1 — A(f>) > Acj> which implies that after a finite number of 
iterations of the return map the firing of the first oscillator does not trigger the avalanche any more. Setting 
C = U (1 - Ac/)) + (<zi — 1)e — 1 condition l|3.34|) is equivalent to 

(3.35) If- 1 (R ((ai - l)e)) - U- 1 (R (()) > U' 1 (1 - (N - a)e) - ET 1 (1 - (N - l)e + Q 

for all (ai — 2)e < ( < (ai — l)e. For £ = (ai — l)e both sides in (|3.35|) are equal and the condition (|3-32|) thus 
ensures ([335]) to hold for all (oi - 2)e < C < (oi - l)e. □ 

Note that for convex U and neuronal R the right hand side of inequality ([3.32P is smaller than one and thus a 
sufficient condition for a\ -avalanches to split up after a finite number of applications of the return map is 

R'(C) > 1 for all C G [(ai - 2)e, (ai - l)e] 

In particular for a partial reset i? c =i(C) = C an avalanches become unstable under the return map and thus by 
corollary [9] only the asynchronous state remains stable for all convex dcpd U and c = 1. 

We used theorem [8] and proposition [15] to determine for a convex LIF rise function I^lif ( c ^ IB.5I eq. (|B.7P ) and 
linear partial reset R c the regime where avalanches of different sizes become unstable under return. The most strict 
condition in (|3.29j) is for a — 1 which yields an implicit equation for the lower bounds on the critical c values below 
which the invariance of di-avalanches is ensured. The upper bound is obtained by (|3.30p using a = 1 and is very 
close to the bound given in proposition [HI The bounds are plotted in fig. 13.71 and are in good agreement with the 
numerical data. 

Near the lower transition point c^ t the system shows aperiodic behavior when starting close to the synchronous 
state. A possible explanation for this dynamics is the competition of two counteracting mechanisms: (i) Large 
avalanches become unstable under return and thus tend to desynchronize the phases which results in a split of the 
avalanche into smaller stable avalanches, (ii) The solution to equation (|3.12p for these asynchronously firing smaller 
clusters involves a* < 0, i.e. the smaller avalanches tend to absorb each other and resynchronize the system yielding 
again larger unstable avalanches. Note that here irregular dynamics arise via a mechanism different form as for 
example network heterogeneity [15] or using excitatory and inhibitory interactions |67| . 

3.6. Extensive Sequence of Desynchronizing Bifurcations A Solvable Example. Figure f3T6l (d-f ) illus- 
trates that the rise function J7& is both icpd and dcpd. In fact, 

AH b ((f), A<f>, e) := H b (0 + A0, e) - H b (0, e) - A(f,e bE 

is independent of <j> and hence 

-^-AH b {4>,Acp,e) =0. 



a a b 




Figure 3.7. Sequential desynchronization in a network (N = 100) with icpd rise function 
(E eq = 1.1, E syn = 3) and linear partial reset R c . (a) Observed cluster sizes of periodic states 
after a time t = 10000. For each c value 100 simulations were started in the synchronous state 
with a small perturbation added. The upper line shows the bounds on a obtained from (|3.30[) in 
theorem [14] above which a-clusters are unstable. The lower line is the bound obtained via i|3.29p 
below which a-clusters are ensured to be stable. The shaded area marks the transition region where 
states other than the synchronous an asynchronous state are observed. In the blue region we find 
no periodic asymptotic dynamics. The dashed lines show the theoretical bounds for the transition 
region, (b) Aperiodic dynamics for c\ = 0.18. 



Thus for Ub equality holds in (|3.26|) and the "best"- and "worst-case" return maps become identical. This property 
allows to obtain exact analytical results. 



Proposition 16. Consider a homogenous excitatory all-to-all network of N pulse-coupled oscillators evolving ac- 
cording to H2.17p - H2.18p with convex rise function Ub (b < 0) and neuronal partial reset Re- 
Then for each 2 < a < N there exist a critical reset strength ci^ such that for all c > c c °^ avalanches of size 
greater or equal to a are unstable under return and avalanches of size smaller than a are stable. For c < all 
avalanches are stable under return. The critical reset strengths are determined from the equation 



( e -b<& 

(3.36) e b(l-[(N-a)+4f(a-l)] e ) = V 



(e~ te - 1) 

and satisfy < c£r < Ccr ^ < • ■ ■ < Ccr < 1. 

Proof. Since Ub is icpd and depd, equality holds in (|3.33p . i.e. for T G T ai 

(3.37) AMjr (A<f>) := 1 — M? (1 — A<b) = 1 - S ai o H {N _ ai)e o J aiE (1 - A0) 

Thus the return map for the phase differences only depends on the avalanche size a\ and is independent of the 
precise form of the other avalanches a, , i > 1 and intermediate shifts <7j . Explicitly 

6e(iV-ai+c(oi-l)) , N 

AMjr (A0) = [e' bc (e» + (l - e b ) A0) c - l) 

for all T G T ai . A straight forward calculation shows that AM? has the properties 

d d 2 

(3.38) AMjr (0) - , — AM r (A0) > and — ^ AM r (A0) < 

Thus if the condition 

(3.39) AMjr (1 - U- 1 (1 - e)) < 1 - U' 1 (1 - e) 

is met all other conditions for 1 < a < ai in p.29p are also satisfied. On the other hand almost all perturbations will 
cause the avalanche to be triggered by a single oscillator. Thus if condition (|3.39[) is not satisfied, i.e. AMj? (A</>) > 
A<p for all A(j) > U^ 1 (1 — e) — 1 the avalanche will split up after a finite number of iterations of the return map. 
Thus H3.39P is a necessary and sufficient condition for stability of an a-cluster under the return map. We are 
interested in the critical strengths for which an a-cluster becomes unstable and hence we use equality in (|3.39|) 
and basic algebra to obtain the implicit expressions (|3.36p for the . 

Since we have assumed {N — 1) e < 1, & < and c e [0, 1] we see that the left hand side of (|3-36[) lies in the 
interval (0, 1) and decreases monotonically with increasing c. The right hand side is for c = and increases 
monotonically with c until it becomes 1 for c = 1. Thus by continuity for all 2 < a < N there always exist a 
solution < Cc?^ < 1 to this equation. Note that the special case a = 2 is explicitly solvable for and yields 

(3.40) eg) = 1 log (l + e -K^-2) E+6 ^ _ e -te^ 

For fixed < c < 1 the left hand side of <|3.36|) is strict monotonically decreasing as a increases whereas the left 
hand side is independent of a, thus < < Cc^~^ < • • • < Ccr < 1- D 

The theoretical prediction (|3.36|) for the desynchronization transition is plotted in Fig. 13.31 and is in excellent 
agreement with the numerically observed transition. 

Remark 17. Note that (|3.36|) involves all relevant network parameter. In particular, by choosing b — > — oo equation 
(|3.40p shows that c a can be made arbitrarily small. This implies that the entire sequence of desynchronizing 
bifurcations may occur over an arbitrary small interval 



r W r ( 2 ) 



Remark 18. We also remark that the number of bifurcation points in this sequence is N — 1. At each bifurcation 
point Ccr' all periodic states with at least one cluster of size a and all other cluster sizes less or equal to a, i.e. an 
extensive combinatorial number of states, becomes unstable simultaneously. 




Figure 3.8. Synchronization and Desynchronization of avalanches in networks with convex rise 
function and partial reset, (a) Sub-threshold inputs synchronize the oscillators. The phase differ- 
ence of a cluster before pulse reception Acj) + is decreased to A0~ afterwards, i.e. A(j> + < A(f>~ . 
(a) Weak partial reset (e.g. c w for R c ) synchronize phase differences: A0 + < A0~ . (b) Due 
to the convexity of the rise function a strong partial reset (c w 1) expands the phase differences 
Aip + > A(ft~ . Clusters lose stability if mechanism in (c) becomes dominant over the stabilizing 
effect (a). 



The mechanism underlying the desynchronization transition are opposing synchronization and desynchronization 
dynamics in the network as illustrated in Fig. 13.81 Due to the convexity of the rise function (a) sub-threshold inputs 
are always synchronizing and stabilize the avalanche, whereas depending on the strength of the partial reset supra- 
threshold inputs in an avalanche can either (b) synchronize or (c) desynchronize the phases. Thus for a weak 
partial reset (e.g. R c with c s» 0) states with large avalanches are stable. When the partial reset is stronger it 
desynchronizes the cluster and, depending on the avalanche size, it may outweigh the synchronization effect due to 
sub-threshold inputs. Larger avalanches receive less synchronizing sub-threshold input from other oscillators and 
simultaneously produce a larger supra-threshold input than smaller ones. Thus they lose invariance under return 
first when increasing the partial reset strength. 

4. Robustness of the Desynchronization Transition 

The desynchronization transition is robust against structural perturbations in the coupling matrix and the rise 
function U. 

4.1. Coupling Strength Inhomogeneity. The desynchronization transition is robust against perturbations in 
the coupling matrix ey. Our numerical experiments show that the transition is also observed when using coupling 
strengths from a uniform distribution on an interval [e m in,£max] for a interval length Ae = e max — £min as large as 
20% of the average coupling strength e — (e m ax — £min) /2 (cf. fig. 14. ip . When Ae becomes larger usually complex 
spike patterns and non-periodic states are observed. 

In inhomogeneous networks sub-threshold inputs of different strengths desynchronize units initially at the same 
phase. Thus coupling inhomogeneity destabilizes clusters of a given size a. In fact, already the lower bound c^ t 
obtained for homogeneous networks via theorem [14] using the coupling strength e over-estimates the stability of the 
clusters as shown in fig. 14.11 The regime where we observe aperiodic dynamics becomes larger in comparison to 
homogeneous networks with the same average coupling strength (e.g. compare fig. 14.11 and fig. I3.7j) , This is due to 
cluster states in the homogeneous network with asymptotic phases of the clusters which are close to an absorption 
(i.e. there are a* » for some i). A perturbation in the coupling now enables the absorption and the restless 
competition between desynchronization and synchronization (cf. sec. 13. 5p induces the aperiodic dynamics. 

Another effect of inhomogeneous coupling is that units synchronized to fire together in clusters are not phase 
synchronized (cf. insets in fig. 14. ip in the asymptotic dynamics. Also intra-cluster phases of equally sized clusters 
and in particular of the asynchronous state show irregular spacings. 

4.1.1. Sigmoidal Rise Functions. Typically rise functions in biological or physical systems are neither purely concave 
nor purely convex. In particular intrinsic neuronal dynamics are often best described with a sigmoidal rise function. 
The quadratic-integrate-and-fire or exponential-integrate-and-fire neuron [2Ql [23] (cf. also Appendix [Bj constitute 
major examples. In networks with sigmoidal rise functions a combination of the effects inherent to concave and 
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Figure 4.1. Sequential desynchronization transition in an inhomogeneous network. Shown are the 
dynamics of a network of N — 50 units with convex rise-function (E cq = 1.1 and -E syn = 3) 
and linear partial reset R c . (d) coupling strengths are drawn from a random uniform distribution 
from Eij £ [0.009, 0.011] excluding self-coupling (en — 0). Simulations are initialized in a perturbed 
synchronous state, (a) Synchronous firing for c\ = 0.12. Inset: Phase are not synchronized 
due to the coupling inhomogeneity. (b) For C2 = 0.19 the synchronous state becomes unstable 
and smaller avalanches are observed in the asymptotic dynamics Inset: As in the synchronous 
state phases of units in a single avalanche are not phase synchronized, (e) For C3 = 0.16 we do 
observe aperiodic firing until t = 15000. (f) For C4 = 0.4 asynchronous firing is observed in the 
asymptotic dynamics. Inset: Spacings of the phases are irregular due to the network heterogeneity, 
(c) Asymptotic cluster sizes of periodic states for different c values starting from 100 different 
perturbed synchronous states. The lines show the lower and upper estimates for the transition 
obtained from theorem [14] assuming a homogeneous network with coupling strength e = 0.01. The 
shaded blue area indicates that states with aperiodic dynamics until t = 15000 are observed. 



convex rise functions influences the network dynamics: Synchronization of units to larger clusters due to the concave 
part (cf. {HI 36]) and stabilization of states with asynchronously firing clusters due to the convex part (cf. theorem 
[3]) . Using strictly neuronal partial resets numerical studies show that for rise functions with dominant concave part 
synchronized firing of oscillators in the asymptotic state is typically found. In contrast, if the convex part is larger 
it is more likely to find clusters of smaller sizes and the asynchronous state. Indeed, for general rise functions U we 
still obtain the stability matrix A in Q3.16P but the non-zero entries l|3.17p can become larger than 1 in the regime 
where U is concave. Thus if the concave part becomes dominant the eigenvalues are no longer bounded by 1 and 
asynchronous clusters states become unstable. 

In fig. 14. 21 a desynchronization transition for the sigmoidal rise function Uq^ and linear partial reset R c is shown. 
In the synchronous state oscillators do not receive any intermediate sub-threshold pulses between successive firing 




Figure 4.2. Sequential desynchronization transition in networks of neural oscillators with a sig- 
moidal rise function. Shown are the dynamics of a homogeneous network (N — 100, e = 0.002) 
with linear partial reset R c and (e) sigmoidal rise function (E syn = 2, a = — 1, /3 = 1). 

Starting with synchrony and inducing a small perturbation (arrow) the network shows (a) aperi- 
odic dynamics for c = c\ = 0.45, (b) clustering for c = c 2 = 0.46 and (c) asynchronous dynamics 
for c = C3 = 0.54. Note the oscillations of the phase which do not appear for purely convex rise 
functions (cf. fig. 13. 7p . (d) Cluster sizes of periodic states observed in the dynamics at t = 5000 
starting from 200 perturbed synchronous states for each value of c. Shaded area marks the transi- 
tion region with states other than solely synchronous or asynchronous. The blue shaded are marks 
occurrence of aperiodic dynamics. The dashed line indicates the critical c^ft determined from 
(|3.32p for a± — N above which synchronous firing becomes unstable. 



and the return map for an oscillator with phase <j> can be written as 

M {{1 ,..., NM (<t>) = U- 1 (R (U(<f>) + (N- l)e - 1)) + 1 - U- 1 (R ((N - l)e)) 

for any partial reset R and any rise function U. After a perturbation the avalanche is typically triggered by a single 
unit and thus the synchronous state becomes unstable if M^^...^}^}^) < 4> f° r an 4> € [l — U^ 1 (1 — e) , l] which 
yields the condition (|3.32j) for ai = N. This can be used to determine the onset of a desynchronization transition in 
the general case as shown in fig. 14.21 (dashed line). The stability of smaller avalanches a\ < N can still be estimated 
with the help of theorem [14] if the rise function is dcpd but not necessarily convex. Conditions for the sigmoidal 
rise functions Uqw and J7q^ to be dcpd are given in appendix iBl 

Desynchronization due to a partial reset has three components: Translation of phase differences into potential 
differences via the rise function U, the relative change of potential differences due to the partial reset R after supra- 
threshold excitation and back-translation of this potential difference into phase differences via C/ _1 (cf. fig. 13.8b ). 
For convex rise functions the slope in the reset zone I R = [0,C/ _1 (R((N — l)e))] is always smaller than in the 
supra-threshold zone I T = [U~ x (1 — (N — l)e) , l] . As a consequence the phase differences in I T are translated via 
U to larger potential differences and the potential differences after reset become larger phase differences during the 
back translation U~ x . This causes an effective phase desynchronization even for partial resets that are non-expansive 
as depicted in fig. 13.8b . 

For general rise functions and non-expansive partial resets the destabilization of a cluster state due to a partial 
reset thus can only occur if the slopes in I T are sufficiently larger than in I R . In fact if this ratio becomes to small 
the transition may not be observed completely for non-expansive partial reset, e.g. for R c in the range c G [0, 1] 
and can be shifted to partial resets that have to be expansive (e.g. for c > 1) . 

Finally note that, in contrast to convex rise functions, for sigmoidal rise functions we always observe "damped 
oscillations" in the Poincare phase plots fig. 14.21 (b.c). The amplitudes of these oscillations become larger when the 



slope of the rise function at the point of inflection becomes smaller. We therefore attribute these oscillations to 
sub-threshold inputs received by oscillators near the inflection point of the rise function. 



5. Discussion 

In summary, we proposed a model of pulse-coupled threshold units with partial reset. This partial reset, an 
intrinsic response property of the local units, acts as a desynchronization mechanism in the collective network 
dynamics. It causes an extensive sequence of desynchronizing bifurcations of cluster states networks of pulse-coupled 
oscillators with convex rise function. This sequential desynchronization transition is robust against structural 
perturbations in the coupling strength and variations of the local subthreshold dynamics. 

Previous studies have not particularly focused on the collective implications of partial reset or similar graded 
resets. In network models with pulses that are extended in time typically a full conservation of the input is 
considered [ 64, 66, .28]. Models with instantaneous responses to inputs consider fully dissipative reset (R{() = in 
our model) gU [26l [0 [55l ED [60] , fully conservative reset (R(() =Q [3 [10] as well as both extremes [31] without 
discussing particular consequences of the reset mechanism. Here we closed this gap and showed that in fact the 
reset mechanism plays an important role in synchronization processes. 

Partial reset in pulse-coupled oscillators keeps the collective network dynamics analytically tractable and at the 
same time describes additional, physically or biologically relevant dynamical features of local units. In neurons, 
for instance, synaptic inputs are collected in the dendrite and then transmitted to the cell body (soma). At the 
soma the integration of the membrane potential takes place and spikes are generated. Remaining input charges on 
the dendrite not used to trigger a spike at the soma may therefore contribute to the potential after somatic reset 

EE1 EH El- 

Such features are effectively modeled by the simple partial reset introduced here. In particular, spike time response 
curves (that may be obtained for any tonically firing neuron (151I53J32I25]) enc °de the shortening of the inter-spike 
intervals (ISI) following an excitatory input at different phases of the neural oscillation. An excitatory stimulus 
that causes the neuron to spike will maximally shorten the ISI in which the stimulus is applied. Additionally, the 
second ISI that follows is typically affected as well, e.g. due to compartmental effects. Exactly this shortening of the 
second ISI is characterized by appropriately choosing a partial reset function in our simplified system. The details 
of such a description and consequences for networks of more complicated neuron models are studied separately [36] . 
For instance, networks of two-compartment conductance based neurons indeed exhibit similar desynchronization 
transitions when varying the coupling between soma and dendrite [36] which in our simplified model controls the 
partial reset. 

The desynchronization due to the partial reset, i.e. due to local processing of supra-threshold input, differs 
strongly from that induced by previously known mechanisms based on, e.g. heterogeneity, noise, or delayed feedback 
[66, 65, 41, 37, 52, [15]. Possibly, this desynchronization mechanism may also be helpful in modified form to prevent 
synchronization in neural activity such as in Parkinson tremor or in epileptic seizures [59} 158]. 

In this work we developed a partial reset for supra-threshold inputs and considered purely homogeneous and 
globally excitatory coupled systems. For inhibitory couplings one can define a lower threshold [H] below which 
inhibitory inputs becomes less effective, i.e. a partial inhibition. In models of neurons, for instance, this could 
characterize shunting inhibition [3]. If two units simultaneously receive inhibitory inputs below a lower threshold, 
a zero partial inhibition, i.e. setting the state of the units to a fixed lower value, is strongly synchronizing in 
analogy to a full reset after supra-threshold excitation. Our findings suggest that similar to a partial reset a less 
synchronizing non-zero partial inhibition may also have a strong influence on the collective network dynamics. In 
biologically more detailed neuronal network models both excitatory and inhibitory couplings as well as complex 
network topologies play important roles in generating irregular [67] and synchronized spiking dynamics pj. An 
interesting task would therefore be to study the impact of partial resets in such networks. 

This work was supported by the Federal Ministry for Education and Research (BMBF) by a grant number 
019Q0430 to the Bernstein Center for Computational Neuroscience (BCCN) Gottingen and by a grant of the Max 
Planck Society to MT. 



Appendix A. The Enestrom-Kakeya Theorem 

A.l. Spectral- Radius and Matrix-Norm. Let A = ay be a n x n matrix. The spectral radius p of a A is 
defined as [32] 

(A.l) p(A) = max ||Ax|| = max |A,| 

x =1 i=l,.. ,n 



where |||] denotes a norm and {Aj}" =1 are the complex eigenvalues of A. If || — 1| is any matrix norm (see [43]) the 
inequality 

(A.2) p(A) < \\A\\ 

is valid and in fact p(A) = inf \\A\\ where the infimum is taken over all matrix norms [32]. Here we only need the 
maximum-absolute-column-sum norm of A defined as 



(A.3) 



IAII = max 



nax > 

l,...,nf-' 



A.2. Companion Matrices. A (n + 1) x (n + 1) companion matrix C has the standard form 

/ ... -co \ 
1 

(A.4) C = 







-ci 



V 1 -c„ / 

with characteristic polynomial 

(A. 5) Pn+iO) = det (z - C) = d + c x z + ... + c n z n + z n+1 

A.3. The Enestrom-Kakeya Theorem. The Enestrom-Kakeya theorem^ [21 133 123 H 132] 

can be stated in the 

following form 

Theorem 19. Let p n (z) — Sj=o c i 7 -^ w ^ c j • > ^ en / or a ^ ^ w ^ M = 



A < max 



Proof. Note first that (3 > 0. We set 



(A.6) 
where 



0<z<n I Ci-)-l 



(z-l)PnW 



Cj-X—PCj 

-co 
c„/3" 



=:/3 



i=0 



1 < i < n 
i = 



Using the definition of (3 one observes that Cj < 0. Comparing (|A.5|1 with (| A.6[) the companion matrix of p n +i is 
given by (| A.4[) . Since 1 + X}j=i = Pn+l(l) = if follows that ||C|| = Y^j=i = ~ Ej=i = 1 when using the 
maximum-absolute-column-sum norm (| A.3[) and hence from (|A.2|) 

P(C) < 1 



Thus for all A with p n +i ( A ) = we have 



p„+i that p„+i 



for A = | and thus |A| < /?. 

Corollary 20. Lei A be a matrix of the form (cf. (|3- 16[) ) 

( -a n ai 
-a„ a 2 



< < 1. For a A with p n (A) = it follows from the definition of 

□ 



A = 



-a n 
V -a„ 



\ 





a n _i 
/ 



wii/i a.; > then 



p(A) < max {ai}™ =1 



1 In 1893 the Swedish actuary and mathematics historian Gustaf Enestrom published this result of roots of certain polynomials with 
real coefficients in a paper on pension insurance (in Swedish) [24]. This result is now often called the Enestrom-Kakeya theorem, since 
S. Kakeya published a similar result in 1912-1913 |34| . But Kakeya's theorem contained a mistake, which was corrected by A. Hurwitz 
in 1913 \M- 



Proof. By a permutation of rows and columns we can cast A into a matrix B = 6jj with non-zero entries = <Jj, 

i € {l,...,n — 1} and = — a„, i € This does not change the spectral radius. The similarity 

transformation to C = Q _1 BQ with Q = diag (gi, . . . , gjv-i) and 51 = 1, q% = II}=i a j-,i ^ {2, • ■ • , n} also preserves 
the spectral radius and C has the form of a companion matrix 1|A.4|) with c, = n?=i+i a« > 0, i 6 {0, . . . , n — 1}. 

Thus p(A) = p(C) < max <i<„ {^7} = maxi<i<„ {aj □ 

Appendix B. Rise Functions 

B.l. Rise Functions for Integrate-and-Fire Models. In this section we derive the rise functions for single 
variable models of the form 

j f v = F (v) + I^t) 

We distinguish between potential independent inputs I m (t) — P(t) with P(t) = J2 S (* — an d the conductance 
based approach I m (t) = g syn P(t) (E syn — v(t)), E syn > 1. More generally, if I- m (t) = Q( v (t))P(t) and Q(t) > the 
transformation 

1 r [t) 1 f 1 1 

(B-l) u(t) = — / TTzrdw Af = / Trrrdv 

y 1 M J Q(v) J Q(v) 

yields 

i.e. a potential independent input Thus if the rise function U for ij n (t) = P(t) is known the conductance based 
rise function U is calculated with the help of (|2 . 13[) as 



lnfl- KT 1Nl 



(B.3) U CB (0) 

111 ^ 1 — i/ syn ) 

The leaky- integrate-and-fire (LIF) model [39] is given by F(u) — —giu + 7 ext which yields 
(B.4) C/lif (0) = E cq (1 - e -«*W) 

where T LW = - j- In (1 - E cq ) and £ oq = ^ + £] > 1. This yields 



(B.5) t££ (</>) 



\n(l-E- y lU hlF (0)) 



In (1 - Esyi) 

For the quadratic-integrate-and-fire (QIF) model [23] with -F(u) = §2 {E T — it) {E t — u) + I CKt one obtains for 
I sya (t)=P(t) 

a — tan (arctan (a) — <p (arctan (a) — arctan (/?))) 



(B.6) C/qif (0) 



a — (3 



wherea=^±^, /3 = a-f, 7 = J ^ - {Et - E,f > 0. Hence 



(B.7) t/CB (^) 



cb / j \ _ ln(l - E sy jU Q i F {<j))) 
QIF (<W - In(l-i^) 



Note that depending on the IF model and coupling type convex, concave and sigmoidal shapes are possible (cf. 
tab IB.ip . We remark that as E syn ->oowe recover the potential independent model from the conductance based 
version, i.e. U CB — * U and the conditions for the different properties of U CB become the conditions for U in tab. 

EU 

B.2. Icpd and Dcpd Rise Functions. Usually it is difficult to verify the icpd or dcpd property (fl2"ll of a rise 
function. Here we show that it is closely related to the third derivative of U . 

We first note that AH obeys the relations AH O, Q, e) = and AH (0, A(f>, 0) = Acf> and hence ^ AH ((f), A<j>, 0) = 
and 

Thus U is icpd if 

(B.8) _jg^ Aff(&A&g)>0 for all (<£,A<M eD 
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Table B.l. Properties of different rise functions. r\ = E syn (a — (3). 



Using < instead of > yields an analogous condition for dcpd U. By definition of AH eq. (|F3.8[) yields the condition 

93 AH ((f), A<f>, g) = 3^^^; e »; y/ y 5 +A ^ 



d(f>dedA(j> r ' (0 + A<?!.,e)) s 

t7" (0 + A0) U" (H (cj) + A0, e)) U' (0 + A0) 2 U'" (H (0 + A<ft, e)) 
U'(H{cb + Acb,e)f U'(H(^ + Acb,e)f 

> V {<j),A<t>,E) €V. 

Substituting H ((f) + A(f>, e) — * (f> and <f> + A(f> — > i/" one obtains 
V ; V y " f/' (0) 17' (^) 2 

as a non-local sufficient condition for a rise function to be icpd. The condition for dcpd U is given when replacing 
< by >. 

Now note that if (]B.9[) is satisfied locally for <f> = ip the sign of the derivative 

d ( u"{<t>? U" ^)U"{4>)U'{<P) \ ( n u"W 2 u"'M \ 

is determined by U" (<fi) since the term in brackets on the right hand side at <f> = ip is positive using inequality (|B.9[) 
and V > 0. Hence, if U is concave, a sufficient local condition for a rise function to be icpd is 

U"U)<Q and U'" (0) < 2 U " jf\ VO < 4> < 1 

U' (0J 

Conversely a local condition for a convex rise functions to be dcpd is given by 

U"((j>)>0 and U"'{(j)) > 2 TT ,,> VO < </> < 1 . 

W (<p) 

Different properties of commonly used rise functions are summarized in Table IB. II 
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